Artificial Staggered Magnetic Field for Ultracold Atoms in Optical Lattices 
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A time-dependent optical lattice with staggered particle current in the tight-binding regime was 
considered that can be described by a time-independent effective lattice model with an artificial 
staggered magnetic field. The low energy description of a single-component fermion in this lattice 
at half-filling is provided by two copies of ideal two-dimensional massless Dirac fermions. The 
Dirac cones are generally anisotropic and can be tuned by the external staggered flux <j>. For 
bosons, the staggered flux modifies the single-particle spectrum such that in the weak coupling 
limit, depending on the flux 4>, distinct superfluid phases are realized. Their properties are discussed, 
the nature of the phase transitions between them is establised, and Bogoliubov theory is used to 
determine their excitation spectra. Then the generalized superfluid-Mott-insulator transition is 
studied in the presence of the staggered flux and the complete phase diagram is established. Finally, 
the momentum distribution of the distinct superfluid phases is obtained, which provides a clear 
experimental signature of each phase in ballistic expansion experiments. 



I. INTRODUCTION 

The preparation of clean condensed-matter systems is 
typically limited by disorder resulting from inevitable im- 
purities, and relevant physical parameters often cannot 
be controlled to high precision. In contrast, ultracold 
gases confined in optical lattices can be controlled to per- 
fection, which permits stringent confrontations between 
experiments and many-body theory. A prominent ex- 
ample is the superfluid-Mott-insulator transition in the 
Bosc-Hubbard model in two [1] and three dimensions 
[2], where experiments have provided a unique quanti- 
tative test ground for the respective theoretical predic- 
tions [3. 4]. Recently, major efforts have been focused on 
reaching the quantum degenerate regime of the fermionic 
Hubbard model with ultracold atoms [5], with the hope 
to promote our present understanding of strongly corre- 
lated electronic systems (e.g., high-T c superconductors). 

The remarkable versatility of optical potentials should 
allow for the realization of the exotic physics known to 
occur for lattice electrons in strong magnetic fields. Un- 
til recently, the generation of artificial gauge fields for 
neutral atoms has been limited to spinning up the en- 
tire system, thereby mimicking the Lorentz force as ex- 
perienced by a charge particle subjected to a magnetic 
field [6-9]. For such systems, the regime of strong cor- 
relations has been shown to be very rich [10]. Reach- 
ing this regime, however, remains a technical challenge 
due to the requirement of rotation frequencies on the or- 
der of the trapping frequencies. Realizations of artificial 
gauge fields, which do not rely on large-scale rotations, 
have been proposed in a variety of theoretical works [11- 
14]. The recent experimental demonstration of a light- 
induced artificial magnetic field by Lin et al. [15], and 
most recently its application to excite a vortex lattice, 
has been a first step to overcome the limitations imposed 
by schemes based upon large scale rotation [16]. 

In conventional solids, the creation of a magnetic flux 



strength on the order of a flux quantum $o = h/e 
through a plaquette is a yet-unaccomplished challenge, 
which has impeded access to the rich physical regime 
characterized by the famous Hofstadter butterfly single- 
particle spectrum [17]. Experiments have thus been lim- 
ited to artificial superlattices with lattice constants on 
the order of 100 nm, where in fact indications of a fractal 
energy spectrum could be observed [18]. In optical lat- 
tices, it should be possible to achieve artificial magnetic 
fields with a magnetic length comparable to the lattice 
length scale. Nontrivial topological properties, such as 
the fractional quantum Hall effect [19] and the anoma- 
lous quantum Hall effect [20], may thus become accessi- 
ble. Theoretical studies of optical lattices with an artifi- 
cial uniform magnetic field [21] and their generalization 
to a non-Abclian gauge field [22] , where phenomena such 
as the Escher "staircase" [23] and the Hofstadter "moth" 
[24] are predicted, have attracted broad interest because 
experimental realizations with ultracold atoms may be 
possible. 

Finally, because of their slow motional time scale, cold 
atoms in optical lattices are well suited for precise manip- 
ulation of the lattice dynamics by external driving [25- 
29] . The theoretical predictions of coherent control in an 
optical lattice with a time-periodic optical potential using 
Floquct theory [25] were successfully tested in a recent 
experiment [29]. These studies have shown that a tem- 
poral modulation, which acts to shake the lattice, can be 
used to modify the effective tunneling strength and even 
to tune it into the regime of negative values. Driven tun- 
neling has also been studied for cold atoms subjected to 
double well potentials; phenomena predicted long ago, 
such as coherent destruction of tunneling [30] , have been 
recently observed [31]. 

This paper discusses how driven tunneling can be 
used to generate an artificial staggered magnetic field for 
atoms in a two-dimensional square optical lattice. A de- 
tailed description of the new physical phenomena that 
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arise when the lattice is loaded with bosons is presented, 
thus extending a recent publication [27]. In addition, a 
discussion for fermions is included. The paper is orga- 
nized as follows. In Sec. II, two different methods are 
used to show how the time-dependent problem yields 
a time-independent effective lattice model with a stag- 
gered magnetic field. Section III describes a study of 
the effective Hamiltonian when the lattice is loaded with 
single-component fcrmionic atoms. The low-energy ex- 
citations at half-filling are shown to behave like Dirac 
particles. The anisotropic Dirac cones are discussed, the 
slope of which is tunable via the strength of the staggered 
flux. At 7r flux we obtain the 7r-fLux phase [32]. Sec- 
tion IV contains a study of the generalized Bose-Hubbard 
model in the presence of staggered flux at zero temper- 
ature. In the weak coupling limit, depending on the 
flux <j), distinct supcrfluid phases are realized: a homoge- 
neous zero-momentum supcrfluid; a staggered- vortex su- 
pcrfluid, characterized by a vortex-antivortex lattice with 
one vortex per plaquette; and a staggered-sign supcrfluid 
with an order parameter with opposite sign for adjacent 
lattice sites. The nature of the phase transitions between 
the different supcrfluid phases is established via a Hartree 
ansatz, and their excitation spectra are studied using Bo- 
goliubov theory. In Sec. V, the superfluid-Mott-insulator 
transition in the strong coupling regime is determined in 
two different ways. The staggered flux renormalizes the 
phase boundary, and, thus, the generalized phase dia- 
gram is obtained with respect to the chemical potential, 
the onsite interaction, and the strength of the artificial 
magnetic field. In Sec. VI, the distinct momentum distri- 
butions of the different supcrfluid phases are calculated, 
which allow for their discrimination in standard ballis- 
tic expansion experiments. Finally, the paper closes with 
conclusions in Sec. VII. 



II. TIME-DEPENDENT OPTICAL LATTICE 
AND THE GENERALIZED HUBBARD MODEL 

A. Time-dependent Hubbard model 

The starting point is the proposal of Ref. [26], which 
pointed out that a refined modulation technique can be 
employed to induce an orbital current with a d x i_ y i sym- 
metry in a two-dimensional optical lattice. The optical 
potential takes the form V(r, t) = Vo(r)+Vi (r, t), consist- 
ing of a stationary part Vb(r) and a temporal modulation 
V]_(r,t) with 



(a) 



Vb(r) = -V p(r), 
Vi(r,£) = KV (r)cos(2S(r)-Qt), 

where p(r) — sin 2 (27r:r/A) + sin 2 (27ry/A), 

_ 1 f sin(27rx/A) — sin(27ry/A) 1 



S(r) = tan 



sin(27ra;/A) + sin(27ry/A) J ' 



(1) 
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FIG. 1: (Color online) (a) Schematic of the staggered current 
(black arrows) driven by the time-dependent optical lattice, 
which leads to two inequivalent sublattices A and B. (b) di 
and d2 are the unit vectors of the A sublattice with length 
d — A/v2, and e;, I = 1, 2, 3, 4 are the nearest-neighbor vec- 
tors connecting the sublattices. In (c) and (d), two distinct 
plaquette summation conventions are defined, denoted by 
and in the text. In (c) the lattice is composed by trans- 
lating an elementary plaquette (shaded area) by means of the 
primitive vectors di and &2 of the A sublattice. In (d) pla- 
quettes are translated by the vectors 2ei and 2e2. The four 
corners of an elementary plaquette are numbered consecu- 
tively by 1,2,3,4 as shown. 



quency, and k is a parameter that quantifies the admix- 
ture of the temporal modulation term. It has been shown 
in Ref. [26] that this optical potential can be engineered 
in experiments by superimposing two bichromatic opti- 
cal standing waves such that Vq can be varied between 
zero and hundreds of the recoil energy Er = 2ir 2 h 2 /m\ 2 , 
where m denotes the mass of the atoms and k can be 
adjusted within the interval [0, 1]. The stationary com- 
ponent Vb(r) of the optical potential provides a regular 
square lattice potential with spacing A/2, whereas the 
temporal modulation term V\{r,t) induces local rotation 
around each plaquette with opposite directions for neigh- 
boring plaquettes [see Fig. 1(a)]; that is, it drives a stag- 
gered current that possesses (d x i_ y i )-like symmetry. 

Given the time-dependent confining potential V(r,t), 
the Hamiltonian in the second quantized form describing 
the quantum gas in the ultracold regime can be written 
as 

H(t) = J d 2 r^(r)^-^-V 2 + V(r,t)^(r) 

d 2 r-0 t (r)V' + (r)V'(r)V'(r), (3) 



1 4iTa s h 2 
2~ 



m 



A is the wavelength of the laser light, Vq is the mean well 
depth of the square lattice potential, f2 is the rotation frc- 



where ip(r) describes a bosonic (fcrmionic) field obeying 
commutation (anticommutation) relations, m is the mass 
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of the bosons (fermions), and a s is the s-wave scatter- 
ing length. Following Rcf. [4], we write the atomic field 
operator ip( r ) m terms of the Wannier wave functions 
VK 1 ") = J2i w ( r ~ where Ri denotes the potential 

minima of Vo(r), at which the atoms are localized. The 
corresponding annihilation (creation) operator is denoted 
by at (aj). The staggered rotation yields a decomposi- 
tion of the square lattice into two sublattices A and B 
[see Fig. 1(b)]. The Bravais lattice is then given by one 
of the sublattices A or B and the unit cell is spanned by 
the lattice unit vectors [see Fig. 1(b)] 

di=ei+e 4 , d 2 = ei+e 2 , (4) 

with the lattice constant d = \j\pl. The four vectors 
e; (I = 1, 2, 3, 4), connecting an A site to its four nearest 
neighboring B sites, are defined by 

ei = -e 3 = - x, e 2 = -e 4 = - y, (5) 

where x and y are the unit vectors in x and y directions 
shown in Fig. 1. In order to distinguish the two sublat- 
tices in our description, we introduce two sets of annihila- 
tion (creation) operators, ai (aj) and bi (b\), correspond- 
ing to the operation on site i of the sublatticc A and B, 
respectively. By substituting the Wannier expansion into 
Eq. (3), we obtain the well-known Hubbard model with 
additional time-dependent one-body terms: 

4 

H(t) = M^A+e, +H.C.) 

r£A 1=1 

reA®B 

4 

+Xi sin(m) Y, ^(-l)' +1 (4&r +e , + H.c.) 

r£A 1=1 

+X2 COS(fit) Y ( n r - n r+ei J J ( 6 ) 

where n r is the number operator on site r. The first 
two terms describe the well-known Hubbard model with 
the nearest-neighbor hopping energy Jo and the onsitc 
interaction strength U given in terms of microscopic pa- 
rameters in the standard manner, J = — f d 2 rw*(x + 
A/4, y + \/4){-(h 2 /2m)V 2 + V (r)]w(x - X/4,y + A/4) 
and U a = \/8Tr(h 2 a s /ma z ) J d 2 r\w(r)\ 4 . Here, har- 
monic confinement of the atoms in the third direction 
is assumed with a localization radius a z . The modu- 
lation amplitudes are given by \i = K ^ J d 2 rw*(x + 
A/4, y)[sin 2 (27ra;/A)-cos 2 (27ry/A)]w(a;-A/4,2/) and X2 = 
2kV Jd 2 r |w(r)| 2 cos(27ra;/A)cos(27ry/A). Note that in 
Eq. (6) we have made the assumption that a single-band 
description, with all atoms residing in the lowest band, 
is sufficient. This requires the rotor frequency f2 to be 
detuned from intcrband resonance transitions of the sys- 
tem. 



For the calculations that follow, it is convenient to 
rewrite the Hamiltonian (6) as 

H(t) = H + W(t)+H int , 
Ho = -JoT, T= J2 4 b h 

<i,j> 

W(t) = Qh iQt + Qe-* m , Q = X - (xzN + iXiM^j , 

M = E ("l)' + 1 («^r +ei +H.C.), 

reA,l=l-4 

t£A 

H in t = \u Q J2 n r (n r -l), (7) 

r£A®B 

where Hq is the stationary kinetic term and -Hint is the 
two-body onsite interaction. As is evident from the 
anisotropy of W(t), the time-dependent part of the opti- 
cal potential Vi(r,t) renders the two sublattices inequiv- 
alent: an anisotropic (quadrupolc-like) time modulation 
of the nearest- neighbor hopping (A4) and of the local 
chemical potential (J\f) arises with a ir/2 relative tempo- 
ral phase lag, which introduces an alternating rotational 
sense to adjacent plaquettes. 



B. Effective staggered magnetic field 

In the following, we discuss two different approaches to 
obtain an effective time-independent description of the 
time-dependent Hamiltonian (7). We begin with an ex- 
pansion of the time-evolution operator of the one-body 
Hamiltonian i/is(t) = Hq + W(t) in a Dyson series. 
Making use of its temporal periodicity and neglecting 
higher-order many-body terms, we obtain an effective 
time-independent Hamiltonian, which turns out to be the 
conventional Bosc-Hubbard model with the kinetic term 
rcnormalizcd by a gauge field. The same result is ob- 
tained upon replacing the classical harmonic oscillation 
by an auxiliary bosonic quantum field, which is subse- 
quently integrated out. Finally, we discuss the gauge 
structure underlying the effective Hamiltonian. 



1. Dyson series 

The time-evolution operator for the one-body Hamil- 
tonian, H 1B (t) = H + W(t), is U 1B (t) = 

T{exp[{-i/h) J* dt'HiB (*')]} where T tt denotes the 
time-ordering operation. For times t which are multiples 
of the revolution time r = 27r/f2 of the rotor potential 
(i.e., t — tit with some integer n), the corresponding 
Dyson series (up to second order in —i/K) is calculated 
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MS 



U 1B {t) = 1 + 



dsHi B (s) 



1 - -Hot - ([H , Q-Q f } + [Q, Q f ]) t 



Q{> 2B,> t) , 



(8) 



where 0(> 2B, > t) denotes terms with at least two-body 
character scaling with t or higher powers of t. Note that 
each expansion order of the Dyson series proportional 
to (—i/K) n with n > 2 can contribute terms linear in t, 
however, each exhibiting at least n — 1-body character. 
With Q as defined in Eq. (7) and upon neglecting the 
0(> 2B, > t) many-body terms, we find 



U 1B (t) ~ 1 



- H eff t 



H 



eff _ 



-J T- 



HQ 



[T,M] 



2Ml 



[M,Af] (9) 



Note that the periodicity i?ie(^ + t) = Hm(t) implies 
that J7ib("-t) = [fiB (■?")]"• Thus, it suffices to justify ne- 
glecting the many-body terms in the derivation of .ffg 
for a single revolution time r. By using the commuta- 
tor (anticommutator) relations of the atomic operators 
according to their bosonic (fermionic) nature, after some 
algebra we find [T, M] = and 

[M,Af]=2 ]T (-l) l {atb r+ei -H.c.} , (10) 

re.4, 2 = 1-4 



and thus 

rreff _ _ 

n — 



j £ { e ^-i) ! /4 fl t 6r+ei+ H.c.}, (11) 



r£A,l=l-4 



where J = y/jf+Wjf, <p = 4tan _1 (Wo/Jo) and 
Wo = xiX2/ftQ- As shown, within an effective 
time-independent description, the temporal modulation 
W{t) renormalizes the real isotropic hopping amplitudes 
Jo of the conventional Hubbard model by adding an 
anisotropic imaginary contribution. Note that the value 
of the phase <f> in Eq. (11) is limited to the interval 
[— 2tt, 2ir], since only Jo > can be accessed with the 
rotor technique. As discussed below, the effective Hamil- 
tonian (11) mimics the action upon charged particles of a 
staggered magnetic field alternating in sign for adjacent 
plaqucttcs. A similar effective description has also been 
used to derive a uniform artificial magnetic field in an 
optical lattice in Ref. [19]. 



2. Auxiliary field method 



The harmonic time dependence in the temporal mod- 
ulation W(t) = Q^e int + Qe.- int in Hamiltonian (7) sug- 



gests a quantization procedure that allows the time de- 
pendence of the system to be eliminated. This method, 
which is similar to the "adiabatic elimination" of the ex- 
cited state of a two-level atom coupled to an off-resonant 
light field [33], amounts to replacing the classical oscil- 
lation terms e ±lfi * with auxiliary creation (annihilation) 
operators p* (p) , which obey bosonic commutation rela- 
tions: 



e~ iat -> p, 

AUt , -t 



pt. (12) 

The one-body part H is (t) = H + W(t) of the Hamilto- 
nian (7) thus becomes 



if is = -JoT + Qty +PQ. 



(13) 



Note that the replacements are carried out such that the 
resulting Hamiltonian is written in an inherently Her- 
mitian form. The Heisenberg equation for p then reads 
ih-^p = \p, H\b\ = , where [p,p^\ = 1 has been used. 
Assuming that the evolution of the auxiliary field is en- 
tirely determined by external driving, thus neglecting any 
back action of the atoms upon the rotor potential, leads 
us to write j^p = —iVLp and thus p = The latter 

may be reintroduced into the Hamiltonian (13), yielding 



H u 



JoT- 

A. 

2MI 



2hn 1 1 



2m 



M 2 



(14) 



Comparison with Eq. (9) shows that the same one-body 
term is recovered, which gives rise to the staggered mag- 
netic field. The nonlocal two-body terms (proportional 
to Af 2 and A4 2 ) are artifacts resulting from our inappro- 
priate implicit assumption that all atoms interact with 
the same quantized mode, a scenario not met in exper- 
iments, where the rotor potential is essentially classical. 
Alternatively, the time dependence can be eliminated by 
means of a path integral method, where the operators p 
and p* are treated as c-numbers in a coherent-state repre- 
sentation. The Lagrangian associated with the Hamilto- 
nian (13) contains terms up to quadratic order in the aux- 
iliary quantum field and we may thus integrate them out 
exactly to arrive at the same effective Hamiltonian (14) 
[34]. 



3. Staggered flux 

A particularly intuitive illustration of the structure of 
the effective one-body Hamiltonian (11) is obtained if the 
lattice is composed of plaquettes translated by the prim- 
itive vectors di and d2 of the A sublattice [see Fig. 1(c)]: 

H e J f = - J y^e ?;0/4 (a\b 2 + b\az + a\b A + b\ ai ) 







+H.c. 



(15) 
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Indices 1 — 4 indicate the four corners of a plaquctte num- 
bered in clockwise order, starting with the lower left cor- 
ner. This representation immediately points out that a 
particle hopping around an elementary plaquette picks 
up an Aharonov-Bohm phase <f> with a sign alternating 
across adjacent plaquettes, which is equivalent to the 
presence of a staggered flux with strength <f> (in units 
of the fundamental flux quantum) in each plaquctte. For 
<f> = ±2ir we have one flux quantum per plaquctte. No- 
tice that to realize this situation for condensed matter, 
lattice electrons would require unrealistically large mag- 
netic fields in the 10 2 — 10 3 Tesla range. 

The time-modulation technique used to derive Hamil- 
tonian (15) lets us only access fluxes <j> m the interval 
[— 2ir, 2tt] because Jo > 0. Note, however, that Hamil- 
tonian (15) displays an 87r periodicity with respect to 
</>, which reflects the existence of a second inequivalcnt 
flux domain for <fi 6 [— 4-7T, — 2ir] U [27r,47r]. This domain 
corresponds to negative values of Jo- A Att change of 
(j>, connecting the two domains, reverses the sign of the 
Hamiltonian. The 87r periodicity with respect to <f> is 
in contrast to the case of a uniform magnetic field in 
a lattice, where the flux per plaquette is defined up to 
an integer multiple of 2tt. The staggered flux in general 
breaks time-reversal and inversion symmetries, except for 
the cases of (f> = 2ixn with n € Z, where the hopping am- 
plitudes attain real (n = even) or imaginary (n = odd) 
values. 

The complex hopping amplitudes are gauge-dependent 
parameters, whereas the total flux passing through a 
closed path is gauge-invariant. Recall that an arbitrary 
lattice Hamiltonian 



H = 



E 

<i,j> 



Xij c \ c j 



H.c. 



(16) 



with complex nearest-neighbor hopping amplitudes obey- 
ing Xji = Xij is invariant under the local U(l) gauge 
transformation 



c 2 exp[ 



Xij -> Xijex-vWj - °i% 



(17) 



It is interesting to note that by means of a gauge change 
in Eq. (15), one can obtain a new Hamiltonian with a 
periodicity in the flux <f> reduced to 27r: on the plaquette 
at position din + d.2m, n, m € Z the replacements a v — > 
av e -^(v+2n+2m)/4 are mad6; where y g {^2,3,4} and 

the A and B operators are not explicitly distinguished 
here. The resulting gauge-transformed Hamiltonian 



H, 



-j'^2(a\b 2 + b\a 3 



i\bi 



+H.c. 



(18) 



trivially exhibits 27r periodicity with regard to </>. 

To describe a finite system, one conveniently chooses a 
periodic boundary condition where the opposite sides of 
the y/~N x \f~N square lattice are identified. This choice 



results in a topology of a torus for the system consid- 
ered. If we now compare the total flux gained around a 
non-contractable loop on the torus, the two Hamiltoni- 
ans in Eqs. (15) and (18) give rise to distinct physical 
realizations. Although there is no net global flux (or </>/4 
flux) gained in the Hamiltonian (15) for TV even (odd), a 
nonzero global flux of \N(f)/2 is accumulated in the ei 
direction for Hamiltonian (18). Throughout this paper 
we will work in the original gauge of Hamiltonian (15), 
which gives the desired physical realization. 



III. 2D MASSLESS DIRAC FERMIONS WITH 
ANISOTROPY 

We now consider loading the optical potential with 
single-component fermionic atoms. In this case, s-wave 
scattering of the atoms is absent due to the Pauli prin- 
ciple. Furthermore, for the low temperatures considered 
here, higher angular momentum collision channels are 
negligible. The effective Hamiltonian in Eq. (11) thus 
provides a complete description of the system in the 
tight-binding limit, realizing an ideal lattice Fermi gas 
in the presence of a staggered flux cj>. By Fourier trans- 
forming the operators 



a k e 




k-O+e;) 



(19) 



Hamiltonian (11) is expressed in momentum space by 



m ff = - E 

kelBZ 



e k a iA + H.c, 



(20) 



with 



£k = 4 J[cos(0/4) cos(fcid/2) cos(fc 2 rf/2) 
-isin(0/4) sin(fcid/2)sin(fc 2 rf/2)], 



(21) 



and the lattice momentum summation is restricted to 
the first Brillouin zone (1BZ) with k u = k ■ d„/d G 
[— ir/d, ir/d], v £ {1,2}. The total number of lattice sites 
is N = 2Na = 2Nb with Na and Nb denoting the num- 
ber of A- and B-sites, respectively. Upon performing the 
canonical transformation 

«k = ^=|^,(-a k + /3k), fe k = 4^( ak + ^ k )' ( 22 ) 
V2 kk| V2 

Hamiltonian (20) becomes 

H e Q ff = Y, (-kkl4/3 k + |e k |4ak), (23) 

where the single-particle spectrum is given by 

|e k | = 2J[cos 2 (k + d) +cos 2 (fc"d) 

+2cos(0/2)cos(fc + d)cos(/c~d)] 1/2 (24) 
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FIG. 2: (Color online) Single-particle spectra of the ideal 
lattice fermions subjected to different staggered fluxes: (a) 
<j> = 0, (b) = tt/4, (c) cf> = 7T, (d) (f) = 7tt/4, (e) = 2tt. 



with fc ± ee (ki ± & 2 )/2. The energy spectrum (shown in 
Fig. 2 for different flux values 0) consists of an upper and 
a lower band due to the bipartite lattice structure. The 
operators and /3£ create a quasiparticle in the upper 
band with energy |ek| and in the lower band with energy 
— | ek | , respectively. Note that ek shares the 87r periodicity 
of Hamiltonian (15) with respect to <p, whereas the energy 
spectrum |ek| exhibits a 4ir periodicity. 

For zero flux <p = 0, the upper and lower energy bands 
recombine at the Brillouin zone edges. By mapping the 
upper energy band to the second Brillouin zone, we re- 
cover the standard tight-binding energy dispersion in the 
absence of a gauge field, with the unit cell consisting of 
an elementary plaquette (the new Brillouin zone is then 
rotated by n /4 and expanded by a factor of V2). The 
presence of the staggered flux immediately leads to in- 
teresting properties in the energy band structure. For 
<p 7^ 2mr, n E Z, the upper and lower energy bands in- 
tersect at four conical points (±7r/d, 0) and (0,±7r/d) on 
the Brillouin zone edges. However, there are only two 
incquivalent points, which we denote K+ = (ir/d,0) and 
K- = (0,7r/(i), given by the zeros of the energy spec- 
trum |e_R- ± | = 0. At half-filling, the lower energy band is 
completely filled and the Fermi level coincides with the 
conical points giving rise to exact particle-hole symmetry. 
An expansion of the energy dispersion for small momenta 



around cither of the conical points K± gives 
kif ± +k| - \/2Jd|[l ± cos(0/2)] k\ 

\ 1/2 

+ [I T cos(0/2)] k\\ + C (|k| 2 ) . (25) 

We see that the low-energy excitations disperse linearly 
in momentum, (i.e., they are Dirac-like) in contrast to the 
case of ordinary particles with a quadratic dispersion. By 
defining the Fermi velocity hvp = \/2Jd, the low-energy 
Hamiltonian becomes 

if e// ~ V2Hv F { [cos(0/4)fc! - i sin(<£/4)fc 2 ] a\ M b +M 

+ [cos(0/4)fc 2 - ism(<f>/4)kx) al k &_, k + H.c.j, (26) 

which contains two copies of Dirac-like particles de- 
scribed by the operators (a + .k,fr+.k) and (<J-,k, fr-.k), 
one around each individual Dirac point K± , respectively. 
Notice that several remarkable phenomena, for example, 
the Klein paradox and the phenomenon of Zitterbewe- 
gung, expected for noninteracting Dirac particles in two 
dimensions, are to be met here. We refer the interested 
reader to the review work about graphene, the prototyp- 
ical system exhibiting Dirac electrons, in Ref. [35]. 

The Dirac cones arising here are generally anisotropic 
(cf. Fig. 2), which results in anisotropic propagation ve- 
locities. The anisotropy of the cone is controlled by the 
staggered flux. Only at the special value (j> = tt do the 
Dirac cones become isotropic. At this point, the system 
simulates the mean-field Hamiltonian of the 7r-flux phase 
proposed by Affleck and Marston to describe the pseu- 
dogap regime of the high-T c cuprates [32]. Furthermore, 
the picture also becomes reminiscent of graphene tight- 
binding physics. 

The adjustable anisotropy of the Dirac cones is a spe- 
cific feature of the staggered-flux scenario in a square 
lattice and is intimately connected to the breaking of 
time-reversal and inversion symmetries. It does not arise 
in graphene or graphcnc-like system with cold atoms in 
a hexagonal optical lattice [36]. For graphene, the inser- 
tion of a time-reversal symmetry breaking perturbation 
would move the two inequivalent Dirac points towards 
each other, or produce a gap in the spectrum, while the 
isotropy of the cones is maintained [37]. According to 
Ref. [38]. anisotropic Dirac cones could be engineered in 
graphene by growing the graphene on top of a suitably 
patterned periodic potential. However, the anisotropy 
would then be fixed. The in situ tuning of the cone 
anisotropy in our system is reminiscent of options aris- 
ing in organic compounds (see Ref. [39]), which provide 
Dirac cones with a tunable tilt. 

In summary, by loading the staggered optical lattice 
with single-component fermions, we obtain an ideal Dirac 
system with tunable anisotropic Dirac cones at half- 
filling. The next section discusses the case for bosons, 
where interactions become important. 
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IV. BOSONIC SUPERFLUID STATES 

For single-component bosons, the non-vanishing s- 
wave collisions between the atoms give rise to the onsitc 
Hubbard interaction H mt of Eq. (7). The staggered flux 
modifies the hopping term of the conventional Hubbard 
model according to Eq. (11). Therefore, we now study 
the generalized Bosc-Hubbard model 

H BH = H e ff + n r {n r -l), (27) 

with Hq given by the Hamiltonian (11). This sec- 
tion considers the weakly interacting regime governed by 
the physics of Bose-Einstein condensation. It is shown 
that, for different flux values 0, distinct superfluid phases 
can be realized: a homogeneous zero-momentum super- 
fluid for — tt < (f> < tt; a finite-momentum superfluid for 
— 37r < (j) < — tt or tt < (j) < 37T, characterized by a 
vortex-antivortex lattice with one vortex per plaquette 
and different rotational directions for the two flux inter- 
vals; and finally, for — Air < (f> < — 3ir and 3tt < <fi < Att, 
a finite-momentum superfluid with an order parameter 
that has opposite sign for adjacent lattice sites. 

For sufficiently weak interactions and low tempera- 
tures, the atoms Bose-condense in the lowest-energy 
single-particle state. The many-body ground state is 
then well described by the Hartree expression |^ko) = 
(Pk ) N ° |0), where No is the number of condensed atoms, 
k is the quasimomentum of the lowest-energy single- 
particle state, and /3^ o is the corresponding quasiparticle 
creation operator introduced in Eq. (22). 

As illustrated in Fig. 3, the minima of the lower band 
of the single-particle spectrum — |ek| in Eq. (24) arise at 
positions in fc-space depending on the value of the flux 4>. 
Two distinct cases arise: the lowest-energy state occurs 
at the center of the Brillouin zone ko = (0, 0) = if 
— tt < (j) + 47rmo < 7r, or it occurs at the four corners 
of the Brillouin zone ko = (±ir/d, ±ir/d) = tt if tt < 
(/j + ATrm^ < Sit. Here, mo and m n are arbitrary integers. 
In the k = case, using e k=0 /l e k=o| = s g n [cos(</'/4)] in 
Eq. (22), we may calculate l^ko) m configuration space, 

l*o,(-i)»o) = 




Depending upon whether nio is even or odd, we obtain 
different superfluid phases. For even mo, which corre- 
sponds to positive values of the hopping strength Jo, we 
recover the familiar zero-momentum homogeneous super- 
fluid state known from the conventional Bosc-Hubbard 
model. For odd mo, the boson operators occur with dif- 
ferent signs for the two sublattices A and B; that is, the 
order parameter is constant except for a different sign at 
the A and 2?-sites (referred to as staggered-sign super- 
fluid). 



To understand better how the staggered-sign super- 
fluid arises, note that although the energy band structure 
(in Fig. 2) remains invariant under the transformation 

a k a k , 6 k ->• -6 k , (29) 

the corresponding upper- and lower-band states are inter- 
changed. The staggered-sign superfluid and the uniform 
superfluid arc thus distinct, despite the fact that they 
arise for the same lattice momentum k = 0. 

In the case ko = tt, each of the four equivalent min- 
ima at the corners of the Brillouin zone, which are re- 
lated to each other by reciprocal lattice vectors, yields 
e k=7 r/l £ k=7r| = i sgn[sin(</>/4)]. Introducing this into 
Eq. (22) leads to the fc-space expression 

l^.(-i)-) = |^((-ir-iat+&t)j | >. 

(30) 

After Fourier-transforming the creation operators a\ = 
N A 1/2 E m ^*U n f™ {m+n) and similarly b\ = 

mdi + nd.2, we obtain the ground-state wave function 
(30) in real space 




where 53 rj denotes the summation over the shaded pla- 
quettes shown in Fig. 1(d). One recognizes that this wave 
function (referred to as staggered-vortex superfluid) ac- 
cumulates a phase of ±2tt, when moving around an el- 
ementary plaquette, with alternating sign for adjacent 
plaqucttes. This forms a lattice of singly quantized stag- 
gered vortices, which are commensurate with the exter- 
nal staggered flux. The Bose-Einstein condensate (BEC) 
formed for the magnetic flux tt < 4> + 47rm 7r < 37r is thus 
characterized by a vortex-antivortex lattice, whereas the 
rotational direction on a given plaquette is determined 
upon whether m w is even or odd. 

Next, we consider the stability of the two BECs which 
can arise for flux values in the interval [— 2tt, 2tt], given by 
the states |\S r o,(-i) m o) and |\E' Wj (_i) m „) in Eqs. (28) and 
(31) for both mo and m w even. We employ a variational 
approach for the BEC ground state with the ansatz 

|f, <j) = -L(e" l « /2 cos(a)/3t + e^ 2 S in( CT )4)^|0),(32) 
V-'Vo 

where the two variational parameters f and a are to 
be determined by minimizing the ground-state energy at 
zero temperature. With respect to the Hamiltonian (27), 
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FIG. 3: Contour plots of the single-particle spectra for differ- 
ent staggered fluxes. Dark regions indicate low energy, (a) 
For 4> = niir, n integer, the minimum occurs at k = (0,0). 
(b) For 4> = 7r + n 2tt, n integer, degenerate minima occur at 
(0,0) and (±n/d, ±7r/d). (c) For 4> = 2n + nAn, n integer, 
equivalent minima occur at (±7r/d, ±7r/d). 



the variational ground-state energy is calculated to be 



(H 



BHI 



4.7 



= -47V Jsm| - I - 



4 J 

UN (N Q - 1) _ 



UN (N - 1) 4 
Jj cos 4 (ct) 

U(N -1)\ AT , 



N 



2N 



E 



MF- 



N cos 2 (a) 
(33) 



The first observation is that the £ dependence drops out 
completely in the mean-field energy Emf- This can be 
understood at the variational level from the fact that 
the Hamiltonian is not sensitive to the relative phase dif- 
ference between the condensation points. Next, we see 
that for ^ < 7r (7T < ^ 2-7r) the k = uniform 
BEC with do = (the k = it staggered- vortex BEC with 
(To = 7r/2) is indeed the absolute minimum of the mean- 
field energy. Finally, the stability of the different ground 
states is verified by allowing a small deviation e from the 
condensation point k . The variation in energy is then 
given by 



(H 



BH)cro+e 



Emf + 

N U 
AJN 



Ae 2 N J 



V2 



sin 



+ 0(£ 4 ). 



(34) 



Since the quantity in the bracket is positive definite, we 
conclude that the ground state is stable against small 
variations. 



For <j) = 7t, the mean-field energy exhibits two degen- 
erate minima at the two points <7q = and <7rj = 7r/2 sep- 
arated by an energy barrier ~ UN 2 /AN . The absence of 
a £ dependence in the mean-field energy precludes a co- 
herent superposition state of the two condensation points 
at the flux value <f> — 71 ■ It thus suggests that the two 
superfluid phases are separated by a first-order quantum 
phase transition, where the order parameter changes dis- 
continuously across this point. We remark that the Bosc- 
Hubbard model with <j) = it is equivalent to the fully 
frustrated Josephson junction model [40]. 

Having shown the stability of the distinct BEC ground 
states, we now study their excitation spectrum using Bo- 
goliubov theory. We first write Hamiltonian (27) in the 
grand canonical ensemble by introducing a chemical po- 
tential \i: 



Jjeff 



(jtN 



(35) 



(-|e k | - M)4/3k + (|ek| - M)4"k 



and for the interactions 



Hint — T7 



U 
N 



E 

ki+k 2 
=k 3 +k 4 



a k! a k 2 ak 3 a k4 



(36) 



we perform the canonical transformation (22). By iden- 
tifying the condensation mode — > ^/Nq + Pk a > we 
perform the Bogoliubov approximation while keeping the 
fluctuation modes only up to quadratic order. The chem- 
ical potential is chosen such that the terms which are 
linear in the fluctuation vanish; that is, 



(37) 



where uq = Nq/N is the condensate density. After some 
algebra, we obtain the action for the fluctuations 



5[$, §+] « - ln Q UN - | ]T <4 • G k 1 • $ k , 



(38) 



k,m 



where <& k = (a k , a_ k , /? k , /3-k) are the fluctuation fields 
and the one-particle Green function Gk in the Nambu 
space is given by 



J 



-M^ 1 = 



in C/A k , ko ^o^k.ko ^ 

ihw m + |e k | + M ko \ n oUB^ M 

^n UB kMa -ihuj. m - |e k | + M ko \n U A kMa 

V h n oUBl ko h n ° UA K* ^ m -|e k |+M ko / 



/-ihhj m + |e k | + Mi 



r 



where M ko = |e ko |+n C/, ^4 k , = B^ = l + cxp(-2iv? k ), S k , = = 1 - cxp(-2iv3 k ), and ^ k = arg(e k ). To 
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obtain the excitation spectrum, we go back from the Mat- 
subara frequency to real time ihu m — > hu and find the 
poles of the 4x4 one-particle Green function. This can 
be easily done by determining the eigenfrequency of the 
equation det[— /iG^ 1 ] = 0, which yields 



hu> = \J |e k | 2 + kk | 2 + 2n £/|e ko | ± 2n U\e k \y/ F kMo , 

where F kiko = cos 2 (^ k ) + 2|e |[|e |/2n C/+l]/n J7. Once 
again, the two branches of the excitation spectrum are 
due to the sublattice degrees of freedom. To examine 
the long wavelength modes, we perform a Taylor expan- 
sion around the condensation momentum kg in the lower 
branch and get Ek m w(k — ko), with the speed of sound 



'<p k d 
v = i I J cos ( — 



<f> k a d 

4 J cos 

4 2 



2n Q U 



corresponding to the Goldstone mode of the broken gauge 
symmetry. 



V. SUPERFLUID-MOTT INSULATING 
TRANSITION 

In this section, we determine the complete phase di- 
agram of the generalized Bose-Hubbard model in the 
strong coupling regime at zero temperature. In the ab- 
sence of the external staggered gauge field, the zero- 
temperature phase diagram of the Bose-Hubbard model 
comprises a superfluid (SF) phase and a Mott insulator 
(MI) phase. These phases are separated by a second- 
order phase transition, driven by quantum fluctuations, 
which is controlled by the dimcnsionlcss number U/AJq. 
When crossing the phase boundary into the SF phase, 
the U(l) gauge symmetry is spontaneously broken, thus 
giving rise to an SF-order parameter. In Sec. IV it 
was shown that, in the presence of the staggered flux 
(j>, the broken-symmetry phase consists of distinct SF 
phases. As the interaction strength is increased, we ex- 
pect a SF-MI transition to take place for each of these 
SF phases. We first use Landau's theory of phase transi- 
tions by introducing a plaquette order parameter, which 
takes into account the various SF phases. Within this 
framework, we determine the critical coupling strength 
(U/4J) C , where the SF order is destroyed. Next, we study 
the Mott regime in detail and derive the excitation spec- 
trum using the path integral formalism. In contrast to 
the Landau theory, we introduce a Hubbard- Stratonovich 
field in the Mott regime to characterize the Mott state 
and treat the hopping terms as perturbations. 



The creation (annihilation) operators a\,{a v ) 1 with v — 
1, 2, 3, 4, are labelled according to the four sites of an el- 
ementary plaquette without explicitly distinguishing A 
and B operators. The Hamiltonian then becomes 

H = -Je^ /4 (aja 2 + a\a 3 + a\a 4 + a\a t ) + H.c. 

u 4 1 
+ -^n,(n,-l) , (39) 

v=l > 



where n„ = aj,a„ and is the summation over the 
shaded plaqucttcs, as shown in Fig. 1(c). We antic- 
ipate broken-symmetry SF phases to emerge for weak 
interactions and introduce a plaquette order parameter 
ip = ("01, ^2, "03) V'4) to characterize them. By performing 
a mean-field decoupling in the hopping term 



a\a v , = (?/>* + at - VC)(VV' + a v > - ip„>) 
~ ip*a v i + alip v > - VCVV', 



(40) 



with z/, z/ <E {1,2,3,4}, we find the mean-field Hamilto- 
nian Hq.mf + Hi^mf in the grand canonical ensemble 



H 



U 



/' 



v=X \ 



Hi. 



MF 







+ (e-^/ 4 V>2 + e^l^)a\ 

+ (e lH H'2 + e^ 4 04)4 + H.c. 



We see that Hq,mf is diagonal in the number-state basis. 
This allows us to calculate the ground-state energy E[ip] 
up to the second order with respect to the perturbation 
H ,mf to get 

E[^\ = n(n- l)U-2np (41) 
+ ^CM^(n,C7,A,0)V^+O(0 4 ), 



A. Landau Theory of Phase Transitions 

For convenience, in this subsection we write the Hamil- 
tonian (27) in the plaquette notation of Fig. 1(c). 



where n is the filling fraction and U = U/4J, p, = fj,/4J 
are the dimcnsionlcss interaction strength and chemical 
potential, respectively. The 4x4 Hermitian matrix M u y 
is given by 
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-10/4 



M(n, U, //, 



i0/4 



E^(n,U,fi} 



E^(n,U,p,) cosO/2) 



i<p/4 



\ 



E(°Xn,U,fi) cos(0/2) 

e-^/ 4 EW(n,U,Ji)cos((j)/2) 



E^(n,U,fi) 



;< ^/ 4 E^(n,U, ft) cos(4>/2) 



-20/4 



,(42) 



where 



E^(n,U,Jl) 



U(n — 1) — \x (j,— Un 



(43) 



In the standard Landau theory, the free energy is ex- 
panded with respect to a scalar order parameter and the 
vanishing of the second-order expansion coefficient de- 
termines the second-order phase transition point. In the 
present extension of Landau's theory, second-order phase 
transitions occur at the zero crossings of the eigenval- 
ues of the matrix M(n,U,p,,<j)) in Eq. (41). There are 
four eigenvectors and respective eigenvalues of the ma- 
trix M(n, U, p,, 4>) corresponding to the four possible SF 
phases found in Eqs. (28) and (31), namely the zero- 
momentum homogeneous SF and the staggered-sign SF, 



Vo,± = (1,±1,1,±1), 



(44) 



£0,± 



2cos((/)/4)[£; (0) (?i,C/,/2)cos((/)/4) ± 1], 



and the two staggered-vortex SF order parameters with 
opposite rotational directions, 



(l,±i,-l,T*), (45) 
2sm((t)/4)[E {0) (n,U,p) sin(0/4) ± 1] . 



Zero crossings exist for £o,+ if — tt < <fi + 47rmo < t and 
Too is even, for eo,_ if — tt < 4> + Amno < tt and too is 
odd, for £tt,+ if 7r < + Amn^ < 3ir and m v is even, 
and for e^,- if 7r < <f> + 47tto w < 3ir and is odd. We 
may thus determine the phase boundaries where there 
is a phase transition between the SF and the MI in the 
different regimes of <p as 



Mo,± 



and 



1 

2 

±- 



1 

2 

±- 



U(2n- 1)tcos( | 



(46) 



U T cos 



U{2n - 1) T sin 



=F AnU cos 



(47) 



U =F sin 



=F 4nC/ sin 



In Fig. 4(a), the surfaces bounding the n = 1 and n = 2 
Mott lobes, given by Eqs. (46) and (47), are shown 
as a function of the experimentally relevant parameters 
(U/iJ , /j,/4Jq, Wo/ Jq) for the first quadrant of the com- 
plex (J +iW )-plane (0 < J ,Wo, i.e., < < 2n). 
Outside the Mott lobes, the two types of SF orders are 
separated by the horizontal plane at Wo = Jo, corre- 
sponding to a flux cf> = 7T. The plane spanned by the 
Wo axis and the white dashed line in Fig. 4(a), given by 
jji/U = 2 — \pl, corresponds to a filling factor of unity. 
For this plane, the complete phase diagram covering the 
entire range [— 47r,47r] of 4> is plotted in Fig. 4(b) . 



B. Effective Action for the Mott State 



We now employ a path integral formulation to derive 
the excitation spectrum of the Mott state in the strong 
coupling regime, thus generalizing a method presented 
in Rcf. [41]. We first write the partition function for the 
generalized Bosc-Hubbard model in terms of the path 
integral Z = J Va*Va exp{— S[a*, a]/h}, where the Eu- 
clidean action in the grand canonical ensemble is given 
by 



S[a*,a]= / dr a,*(r)(fi,9 T — /i)a;(r) 

Jo HeA®B 



<i,j> 



i£A®B 



and the hopping matrix elements Xij are given by 
Xr,r±ei = J exp(j</>/4) and Xr,r±e 2 = J exp(— i0/4) where 
r G A. Since we are interested in the Mott regime 
where onsite interactions are important, we seek to treat 
the hopping terms J2<i •;> Xij a *i T ) a j( T ) as perturba- 
tions. This is achieved by introducing the Hubbard- 
Stratonovich field (tf>i(r), i/>*(t)), such that the hopping 
terms can be decoupled in the following way: 
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z 



J Vip*V^Va*Va&q?^~j dr ^ (i>t{r) - a*{r 



S[a*,a] 



Dip Vipexp 
x / Va*Vaexp 



<i,J> 



<i,]> 



exp 



S [a*,a] 



where the local action So [a* , a] is given by 

S [a*,a}= dr ^ a*(r)(/i<9 T - ^)ai{r) 

+ 2 U a i( T ) a i( T ) a i( T ) a i( T )- 

We then make use of the cumulant expansion formula 

to expand the partition function Z = 
JVip*Vipexp{(-l/h)S eff [ip*,ip}} in powers of 
(-0i(r), -0*(t)) to obtain the effective action S e // [i/j* , tp] ■ 
Here, the expectation value of the field (Ai) g Q taken with 
respect to the weight exp{(— l/ft)So[a*, a}} is defined in 
the usual way: 

(Ai) §0 = 



T>a*T>aAi cxp^ ——So [a*, a] 



Close to the phase transition, where the Mott field van- 
ishes, we keep only terms up to quadratic order in the 
cumulant expansion to get 



Seff[1>%1>)* fdrJ2 

J <i,3> 



21 1 



We note that the expectation values with odd numbers 
of fields (a*,a,i) vanish. Furthermore, the local nature of 
the action So results in the identities 

(ai(r)a*(r')) §0 =5 i:i (a(r)a*(r'))s , 
(^(r) % (r'))s = (a*(r)a*(r')} 3o = 0. 



By going to the momentum space, where the Mott fields 
are expressed as ipi^j[( T ) = ^ k cl^t) exp(ik • r^) and 
^ies( T ) = & k(r) expfik ■ (r 4 + ei)]), and simplifying, 
the effective action becomes 



S eff [r,Tp] = ~ [ drY, 

J 1, 



± W drdr'(a(r)a*(r')) Sa 

1_ J 



a* k (r)a k (T>) + b* k (r)& k (r') |e k | 



Now, since the Mott state with vanishing hopping is 
spanned by Fock states with fixed number of particles, 
the two-points Green function (a(r)a* (V)) can be eval- 
uated exactly to yield 

(a(r)a*(r')) So = 9(r r')(n + l) e -(-M+ntO(r-r ')/» 
+6(t' - T )ne(-^ n - 1 M ( - T '- T) / h . 

Substituting this expression into the effective action, ex- 
panding the fields in the Matsubara frequencies 



v r m 



and using a representation for the step function 

"°° ds e - l < T - T "> 



(48) 



9(r-r') = - f 

J — ! 



2iri <; + ii] 



we finally obtain the effective action up to quadratic or- 
der 
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S eff[o-i,u Jm ^^ UJml bl UJm ,b ktUJm ] - -} 

= E 

k,m 



t 



«k,i 
&k.<. 



where 



/i + nC/ ifc m + — (n — 1)U 



In order to determine the excitation spectrum, we per- 
form an analytic continuation in the frequency space 
ihw m — > hu and locate the poles of the Green func- 
tion G(k, iuj m ). In this case, it amounts to solving 
det[G _1 ] = 0, or 



| £k | 2 (|e k | 2 / 2 -l)=0. 



(49) 



We then obtain two branches of the quasiparticle and 
quasihole spectra in the Mott state, 



-|e k |-2 M +(2n-l)E7 



±Vkk| 2 -(4n + 2)|e k |£/ + ^ 



qp.ph 



|e k |-2 M +(2n-l)E/ 



±y/\e k \ 2 + (4n + 2)\e k \U + U 2 



Since the quasiparticle and quasihole are produced pair- 
wise in the Mott state, we look for the difference in the 
quasiparticle-quasihole spectra to obtain the excitation 
spectrum 



e k | 2 -(47i + 2)|e k |?7 + ;7 2 



(50) 



where the single-particle spectrum |e k | depends implicitly 
on the staggered-flux strength. At a fixed filling, the 
SF-MI transition is then located at the point where the 
gap vanishes. Hence, by fixing n = 1 and evaluating 
_E k = we find the boundaries between the SF and the 
MI phases [see Fig. 4(b)]. Thus, the SF-MI transition has 
been generalized to the case where the critical coupling 
(U/iJ) c also depends on the strength of the staggered 
flux. 



VI. EXPERIMENTAL SIGNATURES OF THE 
DISTINCT SUPERFLUIDS 

A simple method to distinguish the bosonic superflu- 
ids experimentally is to image their momentum distri- 
butions. This is achieved by allowing the system to ex- 
pand ballistically after turning off the confining potential 




J <0 


J >0 


J <0 




<|> = 4 arctanlWo/JJ 





FIG. 4: (Color online) (a) Phase diagram of the generalized 
Bose-Hubbard model subjected to a staggered flux <f> for the 
first quadrant of the complex (Jo + iWo)-plane (0 < Jo, Wo, 
i.e., < <j) < 2ir). Outside the Mott lobes, two types of 
superfluid orders arise, separated by the horizontal plane at 
Wo = Jo corresponding to a flux <f> = n. The plane spanned by 
the Wo axis and the white dashed line, given 
corresponds to unity filling factor, (b) Phase diagram for 
unity filling factor covering the entire allowed range [— 4-7T, 47r] 
of <f>. The range cj> £ [ — 27r, 2tt\, corresponding to positive Jo, 
is accessible by the time-modulation technique discussed in 
Sec. II, which yields a flux <f> = 4arctan[14 7 o/Jo]- The super- 
fluid phases are indicated by their plaquette order parameters 
according to Eqs. (44) and (45). 



and subsequently imaging the atomic density with stan- 
dard techniques. For sufficiently long expansion times, 
the atomic density reflects the initial momentum distri- 
bution. The momentum distribution of the condensed 
atoms is given by the quantity 



(* f (k)*(k)> = Kk)| : 



ikR 



v(k). (51) 



The first factor accounts for the Fourier transform of the 
Wannicr function w(k). The second factor is the struc- 
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(a) 



(b) 




FIG. 5: (Color online) Momentum spectra for (a) the uniform 
(1, 1, 1, 1) superfluid, (b) the staggered-vortex (1, ±i, —1, q=z) 
superfluid, (c) and the staggered-sign (1,-1,1,-1) superfluid. 



ture factor of the Bravais lattice spanned by the vectors 
2ei, 2e2; that is, the sum extends over all shaded plaquc- 
ttes according to Fig. 1(d). The third factor, w(k), is the 
form factor of the elementary plaquette defined by 



t/(k) = 



J2 ^-^(ala p 



(52) 



with V v {y = 1, 2, 3, 4) indicating the positions of the four 
sites in the plaquette. The expectation values (a^a M ) can 
be evaluated for either of the wave functions in Eqs. (28) 
and (31). This task is considerably simplified by observ- 
ing that in the limit of large lattices, these wave functions 
(after some algebra) can be expressed as products of co- 
herent states formed at each lattice site, 



I*) = 



(53) 



rgD iy=l,2,3,4 



where i/j„ denotes the respective order parameter from 

denotes a coherent 



Eq. (44) or Eq. (45) and \yjfie l ~<) 
state at corner v of plaquette r with an average n atoms 
and a phase 7. With the help of Eq. (53), the plaquette 
form factors ^o,±(k) for the homogeneous (1, 1, 1, 1) and 
the staggered-sign (1, —1, 1,-1) superfluids, and ^7r,±(k) 
for the staggered-vortex supcrfluids (1, ±i, — 1, =Fi) are 
evaluated to give 



«o,+ (k) 
«o,-(k) 
?V±(k) 



4cos 2 [ fc r — ] cos 2 ( /c„ — 



4 sin 2 f k x ^- 1 sin 2 ( k * 



47' 



sin 2 f(k x + k y )^j 



+ sin ( (k x - k y )- 



(54) 



with k x = k- i, k y = k- y. The resulting momentum spec- 
tra are shown in Fig. 5. One recognizes the absence of 
the zero momentum peak for the staggered-vortex and 
the staggered-sign supernuids in Figs. 5(b) and 5(c). 
Whereas for the uniform phase lattice momentum and 
momentum are equal, the staggered-sign phase is com- 
posed of momentum components which differ from k = 
by a primitive vector of the reciprocal lattice. The clearly 
different patterns of Bragg peaks permit a direct identi- 
fication of the respective superfluid in experiments. 



VII. DISCUSSIONS AND CONCLUSIONS 

In this paper, a tight-binding model was studied pro- 
vided by an optical square lattice subjected to a time- 
dependent modulation, which excites staggered currents. 
Two different methods were used to show that the time- 
independent effective description of the system is equiva- 
lent to the Hubbard model in the presence of a staggered 
magnetic field. Due to the sublatticc degrees of freedom, 
the single-particle spectrum of the model presents sev- 
eral interesting features, such as two inequivalent conical 
points in the energy band and distinct energy minima 
that depend on the magnitude of the staggered magnetic 
field. Then two cases were considered, first an optical 
lattice loaded with spinless fcrmions and then a lattice 
loaded with bosons. 

When the optical lattice is half-filled with spinless 
fcrmions, the low energy excitations are governed by 
a Dirac-likc dispersion. The problem is then reminis- 
cent of graphenc. However, here the cones are in gen- 
eral anisotropic, with the Fermi velocity controlled by 
the staggered flux. This feature cannot be easily imple- 
mented in graphene, where the Dirac cones arise due to 
the hexagonal lattice geometry. Nevertheless, anisotropic 
Dirac cones can be obtained by growing the graphenc 
layer on top of a periodically patterned potential. This 
method, however, imprints a fixed anisotropy which can- 
not be tuned at will as in the case of the fermionic cold- 
atom system. 

When the optical lattice is loaded with bosons, novel 
superfluid states arise because the location of the mini- 
mum of the single-particle spectrum depends on the stag- 
gered magnetic field: for a staggered flux — tt < <f> < tt, 
the minimum lies at ko = and for weak interactions 
a conventional uniform superfluid phase is realized. For 
tt < <p < 37r, the minimum lies at ko = (±7r/d, ±w/d) 
and thus a finite-momentum superfluid is realized. This 
phase corresponds to a vortex-antivortcx square lattice. 
An analog phase, however, with opposite rotational di- 
rection is realized for — 3tt < <j) < —tt. Finally, for 
— 47T < <f) < —3n and 3ir < <p < 47r, a staggered-sign 
superfluid phase emerges with an order parameter that 
has opposite signs for the two sublattices. We note that 
the time-modulation technique used here to generate the 
staggered magnetic field does not permit to access flux 
values <j) outside the range [— 27r, 27r], where the conven- 
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tional tunnelling strength Jo is positive. It was then 
shown that the different superfluid states are separated 
from each other by first-order phase boundaries within 
the mean-field analysis. For larger U, a second-order 
phase transition to a Mott-insulator arises, where the 
staggered flux renormalizes the critical coupling. Finally, 
the distinct experimental signatures of the supcrfluids 
that could be observed in standard ballistic expansion 
experiments were discussed. 
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